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Abstract 

^_J ■ In this paper we consider a class of robust multilevel precontioners for the Helmholtz equation 

f^ , with high wave number. The key idea in this work is to use the continuous interior penalty finite 

p^ ■ element methods (CIP-FEM) studied in [iniHT] to construct the stable coarse grid correction 

J, I problems. The multilevel methods, based on GMRES smoothing on coarse grids, are then served 

CIh' as a preconditioner in the outer GMRES iteration. In the one dimensional case, convergence 

■^r I property of the modified multilevel methods is analyzed by the local Fourier analysis. From 

our numerical results, we find that the proposed methods are efficient for a reasonable range 
p^ I of frequencies. The performance of the algorithms depends relatively mildly on wave number. 

In particular, only one GMRES smoothing step may guarantee the optimal convergence of our 
^ ^ multilevel algorithm, which remedies the shortcoming of the multilevel algorithm in [5] . 

^^ ' Keywords. Multilevel method, Helmholtz equation, high wave number, continuous penalty finite 

(-H ■ element method, GMRES method, local Fourier analysis 
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1 Introduction 



The efficient and accurate numerical approximation of high frequency wave propagation is of fun- 
!>• ■ damental importance in many applications such as acoustic, electromagnetic, elasticity and geo- 

physical surveys. When the problem is linear and time-harmonic, it can be typically modeled by 
QQ ' Helmholtz equation. The interest of this paper is to consider the Helmholtz equation with Robin 

^«0 ■ boundary condition which is the first order approximation of the radiation condition: 

O' -Au-K^u = f inS7, (1.1) 



dn 



— — \-\Ku = g on dVt, (1-2) 



KA I where VL C M , d = 2, 3 is a polygonal/polyhedral domain, k > is known as the wave number, 

H ' i = \/— T denotes the imaginary unit, and n denotes the unit outward normal to dVL. 

- ■ ■ For high wave number k, the linear system from the discretization of Helmholtz equation is usu- 

ally strongly indefinite, causing most of iterative methods to converge slowly or diverge. In recent 
years, there have been many advances in the development of iterative methods and preconditioners 
for the solution of the Helmholtz equation (cf. [lOpilj). 

Due to high efficiency of multilevel methods for positive definite problems, more and more at- 
tentions have been received to develop robust multilevel methods for Helmholtz equation. However, 
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a direct application of multilevel methods with standard smoothing and coarse grid correction are 
ineffective. Some strategies have been proposed to remedy the problem (cf. [3l[5t ll51 [T6]). For in- 
stance, Elman, Ernst and O'Leary [5| proposed GMRES smoothing together with flexible GMRES 
acceleration. But in order to achieve convergence, a relatively large number of GMRES smoothing 
steps are needed on the intermediate grids. Another approach, the so-called wave-ray multigrid 
methods pi|16j, was proposed by Brandt and Livshits through defining a meaningful coarse problem 
by augmenting the standard V-cycle with ray grids and using coarse grid basis derived from plane 
waves. This method performs well with increasing wave number, but it does not easily general- 
ize to unstructured grids and complicated Helmholtz problems. Alternatively, instead of applying 
multigrid iterations directly to the Helmholtz equation, a class of shifted Laplacian precondition- 
ers [8l[9] has recently attracted a lot of attention, which precondition the Helmholtz equation with a 
complex shifted operator and is shown to be an efficient Krylov method preconditioner. We would 
like to mention that Engquist and Ying [6l[7] recently proposed two new types of sweeping precon- 
ditioners for central difference scheme of the Helmholtz equation based on an approximate LDL^ 
factorization by sweeping the domain layer by layer starting from an absorbing layer. Similar to 
the wave-ray multigrid methods, the new preconditioners have a nearly linear computational cost 
and the number of outer iterations is essentially independent of the number of unknowns and the 
wave number when combined with the GMRES solver. 

Our objective is to develop robust multilevel methods for the Helmholtz problem (jl.ip - ()1.2p . 
Although the pollution error is inherent in the standard finite element methods (FEM) to solve 
Helmholtz equation, the FEM can still be used on fine grid whose mesh size is sufficiently small 
to reduce the pollution error. In a recent work [2]p, the pre-asymptotic analyses of both the FEM 
and the continuous interior penalty finite element methods (CIP-FEM) are given. In particular, 

the well-posedness of standard FEM has been proved under the condition of — < Co(t)^+^, where 
h is mesh size, p is the polynomial order of approximation space, and Cq is a constant independent 
of K, h,p. Thus, without this condition, the well-posedness of standard FEM on coarse grids in the 
multilevel method can not be guaranteed. Moreover, oscillations on the scale of the wavelength 
can not be resolved well by standard FEM on the coarse grids. By contrast, the CIP-FEM has 
been proved in [21] that the associated discrete problem is always well-posed without any mesh 
condition. Intrinsically, the main technique in the stable CIP-FEM is to add a complex shift in the 
bilinear form, which is similar to the idea of shifted Laplace preconditioner approach. Comparing 
to adding a shift to the original problem in shifted Laplace operator, the well-posed CIP-FEM 
is consistent with the original equation and only changes the discrete bilinear form. Based on 
these observations, the new approach proposed in this work is to apply the CIP-FEM to construct 
the stable coarse grid correction problems. Standard Jacobi or Gauss-Seidel smoothers become 
unstable on the coarse grids, especially on the intermediate grids in multilevel methods. Motivated 
by the smoothing approach presented in |5], the smoothing in this work is to use GMRES method 
based on CIP-FEM on the coarse grids, and standard Jacobi or Gauss-Seidel relaxation on the fine 
enough grids. From our numerical experiments, we find that the number of GMRES smoothing 
steps in our algorithm can be much smaller than that in [5J , even if one GMRES smoothing step 
may guarantee the optimal convergence of our multilevel algorithm. 

Our main tool to analyze the multilevel methods for Helmholtz equation is the Local Fourier 
analysis (LFA), which has been introduced for multigrid analysis by Achi Brandt in 1977 ^. 
Comprehensive surveys can be found in |18j and the references therein. We mainly utilize the 
LFA to analyze smoothing properties of relaxations and convergence properties of two- and three- 
level methods in one dimensional case. This may provide quantitative insights into the multilevel 
methods for Helmholtz problem (|l.l|) - (|1.2p . 

The remaining part of this paper is organized as follows: In section [2l we introduce some 



notation, recall the formulation of CIP-FEM, and present the multilevel method for the linear 
system from CIP-FEM approximation. Section [3] is to present the modified multilevel method for 
Helmholtz problem. Standard FEM is used on fine enough grids, on coarse grids the CIP-FEM 
is utilized instead. Section H] is devoted to the LEA of the multilevel method for one dimensional 
Helmholtz problem, we focus on the smoothing analysis and two- and three-level analysis. In the 
last section, we give some numerical results to illustrate the performance of the proposed multilevel 
methods. 

2 Notations and Preliminaries 

2.1 Formulation of CIP-FEM 

Let 7/t be a conforming quasi- uniform triangulation of il, and denote the collection of edges/faces by 
£h, while the set of interior edges/faces by £[ and the set of boundary edges/faces by Sj^ . For any 
T £ Th, we define hx '■= diam(T). Similarly, for e G £hi define he := diam(e). Let h := max^gj-^ hx- 
Throughout this paper we use the standard notations and definitions for Sobolev spaces (cf. [1]). 
In particular, (•, •)q and (•, ■)-^ for E C dQ denote the L^-inner product on complex-valued L'^{Q) 
and i^(5]) spaces, respectively. Denote by (•,•) := {■,-)n and (•,•) := {■,-)dn- 

Now we define the energy space V := H^{il.) n IlrGT H'^i^)- For any v £ V and an interior 
edge/face e = TinT2, where Ti and T2 are two distinct elements of Th with respective outer normals 
ni and n2, we introduce the jump [Vf • n]]|e = 'Vv ■ uiIti + Vu • n2\T2- Define the sesquilinear form 
6/i(-, •) on y X y as follows: 

bfi{u, v) := (Vu, Vv) + J{u, v), u,v G V, 

where 



J(u, v) := Y, nehe {{Vu ■ nl {Vv ■ nj)^ , (2.1) 



where i'je is a complex number with positive imaginary part. The terms in J{u, v) are so-called 
penalty terms and ije are penalty parameters (cf. |19ll21j ). 

Clearly, J{u,v) = if u G H'^{n) and v eV. Thus, if u G H'^{n) is the solution of (frT]) - (fL2]) . 
then there holds 

ahiu,v) ■.= bh{u,v) - k'^{u,v) +in{u,v) = {f,v) + {g,v), v £ V. (2.2) 

We define the CIP approximation space V^ as the standard finite element space of order p, i.e., 

Vh := {vh G H'{n) : v^W G Vp{T), T G %}, 

where Vp{T) is the space of polynomials of degree at most p on T. The CIP finite element approx- 
imation is to find Uh GVh such that 

ahiuh,vh) = {f,Vh) + {g,Vh), Vh£Vh. (2.3) 

It is clear that if the parameters 7e = 0, then the CIP-FEM reduces to the standard FEM. It 
has been proved that the CIP-FEM is stable for any K,h,p > J191I21J . The penalty parameters 
may be tuned to reduce the pollution errors. The numerical results in [19] show that using about 
the same total degrees of freedom (DOFs), the CIP-FEM yields the least phase error comparing to 
the standard FEM and IPDG method [121. 



2.2 Multilevel methods for CIP-FEM 

Let {7I}^Q be a shape regular family of nested geometrically conforming simplicial triangulations 
of the computational domain obtained by successive quasi-uniform refinement of an intentionally 
chosen coarse grid To- We denote by Vi the CIP approximation space on Ti- It is easy to see 
that the spaces {V;}j^q are nested, i.e., Vq C 14 C • • • C Vl- But the bilinear forms {ai{-,-)}i=o 
defined as in ()2.3p on each Vi are nonnested. Thus in this work we consider the following multilevel 
methods for CIP-FEM discretizations, which has also been used for solving the linear systems from 
nonconforming PI finite element approximations (cf. [17j). 

For brevity, we denote by a;(-, •) the bilinear form a/i;(-, •) on V;, where hi is the mesh size of 71. 
Define projections Pf, Qi : Vl ^>- Vi according to 

ai{Pf^'v,w) = aL{v,w), {Qiv,w) = {v,w), Vv eVi, w £Vi. 

The existence and uniqueness of the discrete problem (|2.3p imply the well-posedness of each Pf^ . 
For < / < L, we also define Af : V/ — )• VJ by means of 

{Afv,w) = ai{v,w), yv,wGVi. (2.4) 

Define F; E Vi by 

{Fi,v) = {f,v) + {g,v) yveVi. 

Then the CIP-FEM on level / (cf. ([O])) can be rewritten as: 

Afuf = Fi. (2.5) 

For the smoothing strategy, in fact, when Khi/p is small enough, either of weighted Jacobi and 
Gauss-Seidel relaxation can be applied. Otherwise, GMRES relaxation can be used as a smoother 
and this will be explained in the following sections. To be precise, we describe the smoothing 
strategy as follows: 

Definition 2.1. Let R^ = (^o')"^ Given a > and < i < L, let Sl = {I : Khi/p < a, I < 
I < L} and G^ = {/ : 1 < / < L} \ 5^. If I e Sl, let Rf be the weighted Jacobi relaxation Rf^j or 
Gauss-Seidel relaxation Rf^ based on Af . Otherwise, when / G Gl we use the GMRES relaxation 
based on Ap. 

We remark that we choose a = 0.5 in this paper. This choice is motivated by the efficient 
relaxation of Jacobi and Gauss-Seidel smoothers, and it ensures that the amplification factor for 
the Jacobi and Gauss-Seidel smoothers will not become too large (cf. section 4.3 in the following 
and section 2.1 in ^). When / G Sl, we also define the operator {RfY by 

{{Rffw, v) = [Rfv, w) , yv,we Vl. (2.6) 

We can now state the multilevel method for solving the CIP-FEM discretization system on level 
L which is non-recursive version of multilevel method. 

Algorithm 2.1. Given an arbitrarily chosen initial iterate u^ G Vl, we seek u^ € Vl as follows: 
Let vq = n"~-^. 

1) For / = 0, 1, • • • , L. When I = or I e Sl, 

vi+i = vi + f^iR'rQiiFL - A^vi). 

Otherwise, perform GMRES relaxation for the correction problem Afwi = Qi{Fl — A^vi) 
and set u^+i = vi + HiWi. Here ^i > is a scaling parameter to weaken the influence of the 
error during prolongations. In this paper, we always set m = 0.5. 



2) Forl = L,--- ,1, 0. When I = or I £ Sl, 

V2L+2-1 = V2L+1-1 + MRfYQli^L - ^L^'2L+l-0- 

Otherwise, perform GMRES relaxation based on V2L+1-1 cls in step 1 to get V2l+2-i- 

3) Set n" = V2L+2- 

Here we note that the GMRES relaxation is not hnear in the starting value, and the error 
operator of the above algorithm can not be written directly in a product form which holds only for 
the particular case Gl = 0- When this particular case is considered reasonably, we set 

Ti := fiiRfAfpf and Tf := fniRfYAfPf, I = 0,1, ... ,L, Gl = H), 

Then the error operator of Algorithm 12.11 for the case Gl = can be derived as 

E*mEm, where Em := {I - Tl) ■ ■ ■ {I - Ti){I - To), El, := (/ - ro*)(/ - T^) ...{I-Tl), (2.7) 

where / is the identity operator in Vl- 

3 Modified multilevel methods for standard FEM 

In general, in order to reduce the pollution error, 6-10 grid points per wavelength are typically 
chosen to yield reasonable accuracy. In a unit square domain [0, 1]^, it is well know that nh = 2-k jn^, 
where n^ is the number of points per wavelength (cf. [14J), implies the deterioration of k/i for 
increasing grid points per wavelength. In theory, the well-posedness of discrete solution by CIP- 
FEM holds for any k, h and p, but it is not guaranteed for discrete solution by standard FEM on 
coarser grids with nh/p > G particularly |19ll21j . where the constant G is independent of k, h and 
p. Therefore, in order to establish stable coarse grid correction problems in multilevel methods for 
Helmholtz equation, the CIP-FEM will be applied on the coarse grids. Besides, for standard FEM, 
eigenvalues of discrete system close to the origin may undergo a sign change after discretization 
on a coarser grid. If a sign change occurs, the coarse grid correction does not give a convergence 
acceleration to the finer grid problem but gives a severe convergence degradation instead. This is 
analyzed in [S] and a remedy combining GMRES method as a smoother on coarse grids is proposed. 
This idea has been applied in Algorithm 12.11 and we will also utilize this strategy in the following 
modified multilevel methods to get more efficient smoother for indefinite discrete systems. 

For brevity, let Af denote the linear operator Af^ with the parameters 7e = 0, i.e., the linear 
operator for standard FEM, and energy operator Pf = Pf\^^=Q. The discrete problem on level I 
is Afuf = Fi for CIP-FEM (cf. ([23])) or Afuf = Fi for standard FEM. When the mesh size is 
sufficiently small to reduce the pollution error and satisfy the accuracy requirement, both CIP-FEM 
and standard FEM can be utilized. However, the nonzero elements of linear system from standard 
FEM may be less than that from CIP-FEM, and when the grid is fine enough, the pollution error 
by the standard FEM is also small. Thus, we may still apply the standard FEM on the fine grids. 

Similar to Definition 12. H we use the smoothing strategy on each Vi as follows. 

Definition 3.1. Let R^ = {A^)~^. Given a > and < I < L. lU £ Sl, let Ri be the weighted 
Jacobi relaxation R-^p or Gauss-Seidel relaxation Rf'p based on Af . Otherwise, when I S Gl we 
use the GMRES relaxation based on Af ■ 



When I £ Sl, the operator Rj is defined similarly as in (j2.6p . Now we state the modified 
multilevel method for the discrete system from standard FEM for the Helmholtz problem (jl.ip - 
(fOD as follows. 



Algorithm 3.1. Given an arbitrarily chosen initial iterate u^ £ Vl and integers mi,m2 > 0, we 
seek n" G Vl as follows: 
Let vo = n"~^. 

1) For / = 0, 1, • • • ,L. Given /// > is chosen as in Algorithm ic. 11 when 1 = or I £ Sl, 

vi+i =vi + fiiiRirQiiFL - A^vi), 

where m = mi if I > 0, m = 1 if I = 0. Otherwise, perform mi steps of GMRES relaxation 
for the correction problem Af'wi = Qi{Fl — Af^vi) and set vi^i = vi + iJ-iVJi. 

2) Forl = L,--- ,1, 0. When 1 = orle Sl, 

V2L+2-1 = V2L+1-1 + Hi{RJ)"'Qi{Fl - Af^V2L+l-l), 

where m = nii if I > 0, m = 1 if I = 0. Otherwise, perform m2 steps of GMRES relaxation 
for the correction problem Afwi = Qi{Fl - A^V2l+i-i) and set V2L+2-1 = V2L+1-1 + f^iwi- 

3) Set u"' = V2L+2- 

Remark 3.1. In [5j, in order to prevent unnecessary damping of smoothing modes which should be 
handled by the coarse grid correction, the postsmoothing is always favored over presmoothing (cf. 
section 3.2 in [3]). This is also true in our work. However, due to the utilization of stable coarse 
grid correction method, the above algorithm will converge even when the smoothing is chosen as 
one step in both post- and presmoothing. 

Remark 3.2. Comparing to Algorithm l2.1l the CIP-FEM is only applied in the coarse grid correction 
when / € Gl in the above algorithm. From the first numerical example in this paper, we can see 
that the convergence property of this two algorithms is similar. Actually, this phenomena can also 
be observed in the following LFA. In order to reduce computations, one may prefer Algorithm 13.11 
in practice. 

4 The one dimensional local Fourier analysis 

In this section, we aim to analyze different approaches based on Algorithm 12.11 for the discrete 
system from one dimensional Helmholtz equation, where standard FEM is utilized on the finest 
grid. We focus on the analysis for the discretization from linear continuous interior penalty finite 
element method (CIP-Pl) by the important tool LFA in multigrid analysis. Here we mainly focus 
on the analysis of two-level methods. The LFA of three-level methods is also mentioned. The 
analysis will imply the efficiency of Algorithm 13.11 The following presentation is related to the 
notation and philosophy from |18j . 

4.1 Basic tools in one dimensional local Fourier analysis 

The LFA is based on certain idealized assumptions and simplifications: the boundary conditions are 
neglected and the problem is considered on regular indefinite grids Qh = {x : x = Xj = jh,j S Z}. 
Although the Robin boundary condition ()1.2p and other absorbing boundary conditions are often 

6 



applied in realistic Helmholtz problem, the neglect of boundary condition does usually not affect 
the validity of LFA (cf. [18j). The LFA in this section can be considered as simplification for the 
analysis of multilevel method for Helmholtz problem ()l.mi.2p in one dimension. 

Let Lfi be a discrete operator with a stencil representation L^ = [lj]h,j G ^- For any u^ defined 
on Qh and a fixed point x £ Gh, L^Uh reads in stencil notation as 

LhUh{x) = y^JjUhjx+jh), 

where J C Z is a certain finite index defining the so-called stencil. The basic quantities in the LFA 
are the Fourier modes (ph{6,x) = e ' with x £ Qh and Fourier frequency G M. In fact, the 
frequency 6 can be restricted to the interval {—it, it] as a fact that (ph{0 + 2tt,x) = iph{9,x). It is 
easy to see that the Fourier modes are all the formal eigenfunctions of Lh: 

LhMG, x) = Lh{9)iph{9, x), xegh,Oe (-vr, tt], (4.1) 

where the eigenvalues of Lh can be presented as Lh{9) = Sjgj^jc'^-', which is called the Fourier 
symbol of the operator Lh. Given a so-called low frequency 9^ G 6iow = (— 7r/2,7r/2], its comple- 
mentary frequency 9^ is defined as 

9^ = 9^ - sign(^°)7r. (4.2) 

Interpreting the Fourier modes as coarse grid functions yields 

^h{0\x) = ^2h{29\x) = ^2h{29\x) = 'Ph{9\x), 9^ E Gjow, x € ^ah- 

The Fourier modes iph{9^,x) and iph{9^,x) are called 2/i-harmonics. These Fourier inodes coincide 
on the coarse grid with mesh size H = 2h, and they can be represented by a single coarse grid 
mode ip2h{29^,x). Hence, each low frequency mode is associated with a high frequency mode. For 
a given 9^ G ©lowj define the two dimensional subspace of 2/i-harmonics by 

Efh :=span{<^/,(0O,x),vPft(0\x)}, (4.3) 

where 9^ is defined as in ()4.2p . A crucial observation is that the space £'2/1 is invariant under both 
smoothing operators and correction schemes for general cases by two-level method. The invari- 
ance property holds for many well-known smoothing methods (cf. |18)). such as Jacobi relaxation, 
lexicographical Gauss-Seidel relaxation, et al. 

Let Mh be a discrete two-level operator. In the following, we will show that a block-diagonal 
representation for M^ consists of 2 x 2 blocks Mh{9) (cf. |18)). which denotes the representation 
of Mh on E2h- Then the convergence factor of Mh by the LFA is defined as follows: 

p{Mh) = sxxY>{p{Mh{9)) : 9 G Giow}, 

where p{Mh{9)) is the spectral radius of the matrix Mh{9). We can refer to [IS] for generalizations 
to /c-level analysis. 

4.2 One dimensional Fourier symbols 

Now we give the Fourier symbols of different operators in multilevel method for CIP-Pl discretiza- 
tion of one dimensional Helmholtz equation (jl.ip . We always assume 7e = 7 for some constant 7 
in ()2.ip . Denote by t = nh, R = —1 — 147 ~ ^^/6, 5* = 1 + 187 — t^/3. Since the boundary condition 

7 



is neglected in the LFA, the stencil presentation of discretization operator A^ from (j2.4p can be 
derived as (cf. [20] ) 

A^ = i[i7 R 2S R i^]h. 

Combining the above expression and (j4.ip yields the Fourier symbol of A^ as 

i>) = i(i2,c„s2« + 2iic„s« + 25). (4.4) 

Obviously, the Fourier symbol of standard FEM with piecewise PI approximation (FEM-Pl) is 

For simplicity, we use standard weighted Jacobi (cj-JAC) and lexicographical Gauss-Seidel (GS- 
LEX) relaxations as the smoothers in the LFA. It is easy to derive the weighted Jacobi relaxation 
matrix as S'l^ = I^ — ^D^ A^ , where I^ is indentity matrix, D^ consists of the diagonal of A^ 
and w is a weighted parameter. Due to the fact that D^ = ~^Ih, one can easily deduce the Fourier 
symbol of weighted Jacobi relaxation as follows: 

Si{e) = 1 - '^{i-f cos2e + Rcos9 + S). (4.5) 

The GS-LEX relaxation matrix is Sj^ = {Dh—Lfi)^^Uh, where —Lfi is the strictly lower triangular 
part of Af^ and —Uh is the strictly upper triangular part of A^. The Fourier symbol of Sj^ can 
also be directly derived that 



i?e-'« + i7e-'2e + 25' 



^rw = - .._;i,.L_r2. ,oc - (4-6) 



Note that for the restriction matrix J^ = [r^]^ and x G ^2/1, there holds 

(/2V(r,-))(x) = j;r,ei^-^>,(r,x) = ^r,e'^'''ip2Hi2e',x), a = 0,1. 

By an analogous stencil argument, the stencil presentation of full weighting restriction matrix can 
be derived to be If^ = [1/4, 1/2, 1/4]^'^. Then one can obtain the Fourier symbol of /^ is 

~2h 1 

For the linear prolongation matrix Jg/u one can also derive its Fourier symbol as follows (cf. [18]): 

4.3 Smoothing analysis 

Weighted Jacobi and lexicographical Gauss-Seidel relaxations are the general smoothing operators. 
It is well-known that such two smoothers are unstable especially for linear systems from indefinite 
Helmholtz equation by standard FEM approximations. This is caused by negative eigenvalues 
of the associated linear system and divergence occurs under such two smoothers. To make up 
the problem, an improvement is introduced by the shifted Laplacian preconditioner [HUlOj. which 
preconditions the Helmholtz equation with a complex shifted operator 

-A-(l + i/3)fc2, (4.7) 



where /3 are free parameter. The use of a shift is important to guarantee multilevel convergence, 
whereas the multilevel method for the shifted operator only converges for a sufficiently large shift. 
But this contradicts the fact that the outer Krylov acceleration prefers a small shift. Based on the 
idea of adding a shift to the original problem, we consider the stable CIP-FEM which is consistent 
with original equation and adds a complex shift in the bilinear form. 

Recalling that every two dimensional subspace of 2/i-harmonics £"2^^ with 9^ E ©low is left 
invariant under the w-JAC and GS-LEX relaxations, then the Fourier representation of smoother 
Sh 



:!GS 



?eo 



Sj^ or Sf^^ with respect to £"2^ can be written as 







^ 



(4.8) 



where Sh{0) is the smoother symbol derived in ()4.5|) and ()4.6p . The spectral radius of the smoother 
operator can be easily calculated since the above matrix is diagonal. 

For simplicity, we concern on the weighted Jacobi relaxation. The frequency 9 maximizing 

over (— 7r,7r] can be calculated by its first and second derivatives, which reveal ^ = or 



;-/ 



J^ 



9 = TT maximizing \Si^{9)\. Hence, the spectral radius of S^f^ can be deduced as follows: p{S'^ 
max{|5;,(0)|,|S;,(7r)}|, where 



\siio)\ 



1 



UJ 



"^-Xn + m, \sii7r)\ = \1 - u; - ^[ij 



R)\. 



S' ' '" ' ''' " ' S 

Since the parameter 7 in CIP-FEM may influence the pollution error, it is critical to make a suitable 
choice. For the case t = nh < 1, it has been derived in [20] that for one dimensional problem there 
exists an optimal choice ijo = "'^lofi-* ^°t)^^ * such that the pollution error vanishes when the 



penalty parameter is chosen as |i7 — ijo] ^ 



■zTj-, where the constant C is independent of k and h. 
The left graph of Figure [1] shows the Fourier symbols Sf^ {9) for cj- JAC smoother with 7 = 70 and 

a; = 0.6. We find that Sfi{9) > 1 always occur at the low frequencies, and small t leads to a better 
relaxation. Similar phenomenon is observed for GS-LEX smoother in the right graph of Figure 
[TJ Thus, for fixed wave number k, both li;-JAC and GS-LEX relaxations can be used as smoother 
on fine grid. Actually, the relaxation properties of these two smoothers are similar when 17 is a 
complex number. 





J QQ 

Figure 1: \Sf,{0)\ with co = 0.6 (left) and \S,, {9)\ (right) over (-7r,7r] for t = 0.1,0.5, 1. 



When a;-JAC smoother is utilized in the standard PI FEM approximation, then 



p{S()\j=o = max{ I 1 



uj{l2-t^ 
6-2*2 



1 + 



6-2t2 



}• 



(4.9) 



Thus, for t = nh near vS, the error is ampHfied by the smoothing. For CIP-FEM, there is a 
(complex) shift in the bihnear form, therefore the amphfication factor is reduced. This may permit 
the use of w-JAC smoother again. In fact, in our modified multilevel method, GMRES iteration 
is used on the intermediate grids (cf. [5]). In contrast with a;- J AC and GS-LEX relaxations, the 
Fourier symbol can not be derived for GMRES smoothing. In the following, we will give some 
explanations for the performance of GMRES smoothing from the numerical approach. 




Figure 2: Amplification factor of GS, Jacobi and GMRES smoothing for t — 0.8. 

For simplification, we consider the one dimensional Helmholtz equation on an interval (0, 10) 
with homogeneous Dirichlet boundary conditions. The piecewise PI CIP-FEM discretization leads 
to a linear system with N x N coefficient matrix (cf. [20j ) 



A? 



/ 25 - i7 R i7 

R 2S R ij 

i7 R 2S R i7 



\ 



V 



ij R 2S R i7 
i-f R 2S R 

i7 R 2S -ij J 



NxN 



where A^ = 10/h — 1. For k = 200, we apply the grid with mesh size h = 0.004, i.e., t = 0.8, 
and choose 7 = 70. Let the vector uq = e^^^'^ be an initial choice for smoothing, where x = 
[xi, • • • ,XAr_i],Xfc = Xk-i + h,xo = 0,k = 1,- ■ ■ ,N — 1, 6 £ (— 7r,7r]. We assume that Sh is the 
relaxation iteration matrix and ui = ShUQ is the new vector after one step of smoothing. Then 
for fixed 6, we obtain the amplification factor for one step of smoothing ps{0) = ||^i||/||^o||; where 
II • II stands for the Euclidean norm. Figure [2] shows the amplification factors of three different 
kinds of smoothing strategies: GS, Jacobi and GMRES relaxations. We can see that for high and 
low frequencies 9, the amplification factor of GMRES relaxation is always smaller than that of GS 
and Jacobi relaxations. The GMRES relaxation can still be convergent even when the other two 
smoothers lead to a divergent relaxation. The similar phenomenon can also be observed for other 
choices of t and 7. Therefore, when standard Jacobi or GS relaxation fails as a smoother, we can 
replace this with the GMRES smoothing. 

From the above analysis, we can see that for the relaxation of linear system for Helmholtz 
problem, the standard Jacobi or GS relaxation can be performed on the fine grids, but they fail on 
the coarse grids. In particular, the GMRES relaxation is efficient for smoothing on the intermediate 
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coarse grids. In the next section, the correction scheme will be taken into account to obtain a more 
realistic convergence analysis of the multilevel method. 

4.4 Two- and three-level local Fourier analysis 

For more details about convergence property of multilevel method in Algorithm 13. 11 we will focus on 
the LFA of two-level method. The LFA of three-level method will also be mentioned. For simplicity, 
we consider the case without postsmoothing. Since the Fourier symbol can not be obtained for 
GMRES smoothing, we only consider the two and three level methods with w-JAC or GS-LEX 
relaxation. Then the iteration operator of Algorithm 13.11 in this case is the nonsymmetric version 
and can be derived as in (|2.7|) . Three cases of iteration operator will be analyzed in the following: 
multilevel methods by shifted Laplace approach, multilevel methods with stable CIP-FEM coarse 
grid corrections, and multilevel methods with stable CIP-FEM for fine grid smoothing and coarse 
grid corrections. 

The iteration operator of two-level method can be written as (/ — Ti)(/ — Tq), and the operator 
Ti is chosen for the corresponding algorithm. Since standard FEM is applied on the finest grid, we 
have Ti = fuRiAf Pf , where Ri is smoother determined by the algorithm. When CIP-FEM is used 
for smoothing on the grids 71 and To respectively, the iteration matrix is given by 



M^ 



iC\ 



Mli-S'OiA^) 



-'Af 



f^oiUA^'y 



ItOaF 



Here, for / > 0, 



S? 



R?A? 



is smoothing relaxation matrix, i"; with the same size as 



A; is identity matrix, If{s > I) is prolongation matrix from level / to s, If{s < I) is restriction 
matrix from level / to s, and i?p stands for matrix representation of smoother Rf with CIP-FEM 
approximation. Let M2 (cf- [4J) stand for two- level iteration matrix, taking the smoothing based 
on shifted Laplace operator (|4.7p with standard FEM approximation. Similarly, when applying 
standard FEM for smoothing on 7i and CIP-FEM for correction on To, the associated iteration 
matrix is derived as 



Af^ 



^i(/i-5f))(li-/.o/J(Ao^)-i/?Af 



where S^ = /i — R^ A^ , and R^ stands for smoothing matrix representation of Ri with standard 
FEM approximation. 

Since every two dimensional subspace (|4.3p of 2/i-harmonics -E2/1 with 6^ G (— 7r/2, 7r/2] is left 
invariant under w-JAC or GS-LEX smoothing operator and correction operator, the representation 



of two-level iteration matrix of Af 2 on E'. 



2/11 



is given by a 2 x 2 matrix as follows: 



Mo 



■c 



( 
Il -1^1 

\ 


II- 


's 
s 




dI . 


Af (00) " 
Af (01) _ 


D 


1 


" Af (00) 
. Af (01) 




D 


?i -^0 






Ao{2 


QOyl 


" I°(0O) ■ 


t 




Af (00) " 
Af (01) _ 


D_ 





(4.10) 



where i"i is 2 x 2 identity matrix and the subscript-/) denotes the transformation of a vector 

— SL 

into a diagonal matrix. Similarly, the representation M2 of two-level iteration matrix based on 
shifted Laplace operator, can be easily obtained (cf. [4]). For the iteration operator M2 , its 
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representation on £'2/1 is given by 

7i - ^1 Ii 



Mn 



■FC 



II - fJ'O 






J D/ J 






;C 



Ao (20' 



0\-l 



hie') 



5 (^°) 

Af(0i) 



D 



(4.11) 



— SL — FC — C PI 

Then the spectral radius of M2 , M2 and 7W"2 for different 6 S Oiow can be obtained analyticahy 
and numerically. 




-1.5 -1 




--—' SL 

Figure 3: Left: Case 1-Case 3: ^(Mj ) with /3 = 0.5, 0.3 and 0.1 over 0° G 6iow, Case 4 and Case 5 denote 

— FC — c — c 

for p{M2 ) and p{M2 ) with ij = i7o + 0.01i {ijo « -0.085 when t = 0.8) on the finest grid. Right: piM^ ) 

with different 17. Since t := fcft, > 1 on the coarse grid, we formally choose the corresponding parameter 

17 = i'jo + 0.051 for the cases in this figure. 

1-1 -—SL 

The left graph of Figure [3] shows the spectral radius of M2 with (3 = 0.5,0.3 and 0.1 over 
„ — FC — C 

6 G ©low, and M2 , M2 under w-JAC relaxation for t = 0.8 on the finest grid. In order to 

make some comparisons between different algorithms, the penalty parameter 17 in the CIP-FEM is 

always chosen as a complex number in the following if there is no special notation. If no confusion is 

possible, we will always denote t = constant for t = nhi on the finest grid. The parameter in uj-J AC 

relaxation is always chosen as w = 0.6 in this section. We observe that for most of the frequencies 

9^ G ©low the spectral radius is smaller than one. The spectral radius tends to larger than one only 

under a few frequencies. Actually, the appearance of such a resonance is caused by the coarse grid 



'?■ 



A:{2e<') 



correction and originates from the inversion of the coarse grid discretization symbol, Aq [W) m 

(j4.10p and (j4.1ip for instance. The minimum of \Aq {29^)\ maximizes the spectral radius of M2 

— c 
and M2 for fixed t and 17. 

Moreover, for different choices of parameter /3 in shifted Laplace operator, the properties of 
— FC ~C ... . ~5L 

M2 and M2 with 17 = i'jo + O.Oli on the finest grid always perform better than that of M2 ■ 

Although the spectral radius of M2 with /3 = 0.1 is almost equivalent to that of M2 o^^r the 

frequencies spreading near zero, there is a relatively large resonance causing by the coarse grid 

correction. Small /3 in ()4.7p deteriorates the coarse grid correction. However, large /3 amplifies the 

corresponding spectral radius over the frequencies spreading near zero. We can refer to a recent 

work [4J about the choice of minimal complex shift parameter in shifted Laplace preconditioner. 

— FC — c 
The choice of 17 in M2 or iVf 2 would also influence its spectral radius. The right graph of Figure 
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(J 

[3] shows the spectral radius of M2 with different 17 under w-JAC relaxation for t = 0.8. The 

' — 'C 
influence of 17 in iWg is similar to j3 in shifted Laplace preconditioner. For other small t < 1, the 

— SL — FC -—c 

properties of M.2 , M2 and M2 can also be observed as above. Moreover, when t = Kh on 

— SL — FC — c 

the finest grid is small enough, the spectral radiuses of M2 , M2 and M2 can all be smaller 

than one over all 6^ G ©low Due to the fact that our aim is to study the influences of smoothing 

corrections in two- and three-level methods, in the following we do not always assume t to be 

sufficiently small. 




SL, 



Figure 4: Case 1 (t = ^/3) and Case 3 (i = 4) denote for p{M2 ) with /3 = 0.8 over 6*° e Glow, Case 2 
{t = VS) and Case 4 (i = 4) denote for piM^ ) with 7 = 0.8. 



^ ^ O J , WC - C 

Alternatively, for large i, the maximum of spectral radiuses of M2 ) -^2 ^'^d M2 appear to 
be spread around 9^ = 0. Actually, for large t the main factor determining the spectral radius is 
the smoothing operator rather than the coarse grid correction. We assume that w-JAC relaxation 

is used as smoother. It can be observed from Figure H] that p{M2 ) and p{M2 ) reach a maximum 
in 6**^ = for t = \/3 and t = 4. Indeed, it has been analyzed in previous section and [Ij that p{SJ^) 
for CIP-FEM and shifted Laplace preconditioner (j4.7p with continuous piecewise PI approximation 
reaches a maximum at = or = vr. From ()4.9p . it is easy to see that the a;- J AC relaxation 
for discrete system from standard PI FEM approximation is inefficient when t is near \/3. Thus, 
two-level method M2 is inefficient in this case. However, M2 and iVf 2 can still be applied due 
to the adding of complex parts on the fine and coarse grids. One can observe from Figure H] that 

the spectral radiuses of p{M2 ) ^^^ p{^2 ) i^ always larger than one when t = \/3, which implies 
the divergence of the associated two-level methods when applying w-JAC smoother. However, the 
spectral radius is smaller than one when t = 4, which may permit the use of w- JAC as a smoother 
on very coarse grids, but we shall not use this and utilize GMRES smoothing on coarse grids 
instead. 

Furthermore, to get a comprehensive insight into the influence of multiple coarse grid correc- 
tions, we next carry out a three-level local Fourier analysis. The iteration operator of three-level 
method by nonsymmetric version of Algorithm 13. 1 1 is derived as {I — T2){I — Ti)(I — Tq). Similar to 
the two-level method, when CIP-FEM is used for smoothing on 72, Ti and To, the iteration matrix 
can be deduced to be 



M' 



c 



l2-P2{l2-S^)iA^)-'Ai 



C\ 



iCn-Ii-1 



l2-piIi{Ii-SY)iA^)-'nA'2 



-poiiiA^r'nA^ 



Define the four dimensional 4h-harmonics by 
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where 0"° 



lal 



2 '' 



sign(-2-)7r, a = 0,1. For any 0° S 0iowi the three-level operators leaves 



the space of 4/i-harmonics E^f^ invariant (cf. [18j). This yields a block diagonal representation of 



tC 



iWg with the following 4x4 matrix M^ 



'C 



Mn 



-c 



/^2 U 



Ani"i 



'C 



'-^Z 



{A,) 



'C, 



5f(0i) 



D> 



■C 






:C , 



l2-f^ollllA''o{29^)-Hl')\llYAf 



(4)*4 



D 



(4.12) 



'C 



where 72 is 4x4 identity matrix, ^2 



00 \ 



ij{e 
i{{e 



01 ^ 



10^ 







i?(0iM 



,^0 






^?(0°°) " 




" 1^(000) " 




" If (000) 




,A?^ 


4(0^°) 


,Af^ 


Af(0Oi) 
Af(0iO) 


5?(^") J 


D 


L A? (011) J 


D 


[ A^ (011) 



D 



'1 ~0 



, and 1 2^ 1 1 are defined similarly. For another approach, we 



assume that standard FEM is applied on Ti for smoothing and CIP-FEM is used for correction on 
7i and To, then the associated iteration matrix is given by 



M^ 



^2(^2 



St, 



C\-lrl 



/2 - Milf (Ii - 5^)(A^)-iliA^ /2 - ^lo^{A, 



^2.4C^-lJ0^F 



Moreover, similar to iWf 1 ^^ denote by M^ the three-level iteration matrix, which takes the 
smoothing based on shifted Laplace operator (j4.7p with standard FEM approximation. Then the 

representations M^ , M^ with respect to M^ and iVf 3 respectively on the 4/i-harmonics 
^Ah2 ^^'^ ^^ derived directly. 





over 61° e 9 



low for the cases t 
FC , 



0.6,0.4,0.2. Right (t = 0.4): Case 1-Case 2: 



Figure 5: Left: p{M^ 

—-SL — C 

p{M^ ), p{M^ ), Case 3-Case 4: p{M^ ) with one step and two step intermediate coarse grid (the 2nd 
level) correction. w-JAC relaxation is applied for smoothing. The parameters in shifted Laplace operator 
and CIP-FEM are chosen as /? = 0.5, i'j — ijo + O.Oli when t < 1 and 17 = ijo + 0.05i when t > I. 
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For fixed wave number k, the performance of three-level methods is shown in Figure [5l For 
instance, one can observe from the left graph of Figure [5] that a coarser initial grid used for M^ 
will deteriorate its convergence. The right graph of Figure [5] shows the spectral radiuses for three 
kinds of approaches: M^^, M^ and Mf^. For t = 0.4, both of Mf'^ and Mf perform better 
than M^ over most of the frequencies, and the spectral radius of M^ is similar to that of M'^ 
over most of frequencies in this case. Actually, the performance of iWg is always comparable with 
M'^ . This is also true for two-level method, which can be observed from the left graph of Figure 
[3l Thus, in order to reduce the computational cost, the modified multilevel methods (Algorithm 
13. ip with stable CIP-FEM corrections on coarser grids and standard FEM corrections on finer grids 
can be usually utilized in practical. From the right graph of Figure [5l we can also see that the 
convergence of M^ does not always perform better with more steps of intermediate coarse grid 
correction. 

5 Numerical results 

In this section, we present some numerical results to illustrate the performance of Algorithm 12.11 
and modified multilevel method (Algorithm [3TT]) for two Helmholtz problems in two dimension. The 
multilevel method is used as a preconditioner in outer GMRES iterations (PGMRES). Gauss-Seidel 
relaxation is always used as smoother when I € Sl, and the smoothing steps are always chosen as 
one step if there is no any annotation. At the l-th. level, the discrete problem reads A;U; = F;. We 
denote by r" = F^ — A^u" the residual with respect to the n-th iteration. The PGMRES algorithm 
terminates when 

||rr||/||rO|| < 10^6_ 

The number of iteration steps required to achieve the desired accuracy is denoted by iter. 

Example 5.1. We consider a two dimensional Helmholtz equation with the first order absorbing 
boundary condition (cf. p^[T9]): 

A 2 r sin(Kr) 

— Au — K u = J := m il, 

r 

^" , ■ ^n 

— — h iKU = g on oiZ. 

on 
Here fi is a unit square with center (0, 0) and g is chosen such that the exact solution is 

cos(Kr) cosK-|-isinK 

K k{Jo{k) + iJi{k)) 

where Ju{z) are Bessel functions of the first kind. 

Table 1: Iteration counts of PGMRES for discrete system from CIP-Pl and CIP-P2 (k = 100). The smoothing 
relaxations on fine and coarse grids are all based on CIP-FEM. 



Level 


2 


3 


4 


5 


DOFs 


16641 


66049 


263169 


1050625 


iter (PI) 


27 


24 


21 


20 


iter (P2) 


26 


22 


19 


18 



In this example, we assume that the coarsest level of multilevel method satisfies k/iq/p ~ 2 for 
K < 360. For 400 < k < 600, we choose the coarsest grid condition such that kHq/p ~ 1.1 ~ 1.7. 
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The parameters in dHJ are chosen as 7^ = 0.01 + 0.071 (cf. [I2]) for CIP-Pl, and 7^ = 0.005 + 0.0351 
for piecewise P2 CIP-FEM (CIP-P2). We first test the algorithms for the case with wave number 
K = 100. When Algorithm 12.11 is applied to the discrete system A^ul = Fl, the smoothing 
relaxations are all based on the CIP-FEM approximation on fine and coarse grids. Table [T] shows 
the corresponding iteration counts of PGMRES for discrete system from CIP-Pl and CIP-P2 
approximations. We can observe that for fixed k the algorithm is robust on different levels. 



Table 2: Iteration counts of PGMRES for discrete system from FEM-Pl. The smoothing relaxations on fine and 
coarse grids are all based on FEM-Pl for «: = 100, 200. 



K - 


= 100 


Level 


2 


3 


4 


5 


DOFs 


16641 


66049 


263169 


1050625 


iter (m2 = 1) 


58 


87 


94 


90 


iter (m-z = 20) 


46 


53 


49 


47 


K - 


= 200 


Level 


2 


3 


4 


5 


DOFs 


66049 


263169 


1050625 


4198401 


iter (7n2 = l) 


177 


>200 


>200 


>200 



For fixed k, when the grid is fine enough to get the accuracy, standard FEM can be utilized again 
to discretize the problem. But standard multigrid method fails to solve the corresponding discrete 
system for the case with large wave number. In the following, we mainly apply the multilevel 
method presented in Algorithm 13.11 with different smoothing strategies. Table [5] shows the iteration 
counts of this multilevel-preconditioned GMRES method with GMRES smoothing on coarse grids 
for K = 100, 200, and the smoothing relaxations on fine and coarse grids are all based on standard 
PI FEM (FEM-Pl) approximation. We find that the iteration count is mesh independent for 
fixed K, but it increases rapidly with larger wave number. For instance, for the case k = 200 
with 7712 = 1) the iteration counts of PGMRES will be more than 200, even when the GMRES 
smoothing is performed by m2 = 20 steps, the iteration counts are still more than 100. Although 
more steps of GMRES smoothing can reduce the total iteration counts, it requires more memory 
to store data in the computation. Actually, the slow convergence of the algorithm mainly lies in 
the bad approximation on coarse grids. Next, we apply CIP-FEM to construct stable coarse grid 
corrections. 



Imaginary part of discrete FEM-P2 soiution 



imaginary part of discrete FEM-P2 soiution 





Figure 6: Surface plot of imaginary part of discrete FEM-P2 solutions for n — 100 (left) and k — 200 (right) on the 
grid with mesh condition Kh/p « 0.55. 

Figure [6] displays the surface plot of imaginary part of discrete FEM-P2 solutions for k = 
100, 200 on the grid with mesh condition nh/p ~ 0.55. Indeed, the discrete solutions have correct 
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Table 3: Iteration counts of PGMRES for discrete system from FEM-Pl approximation (k — 100). Case 1: the 
smoothing relaxations on fine and coarse grids are both based on shifted Laplace operator H4.7|) with /3 = 0.2 by 
FEM-Pl, Case 2: the smoothing relaxations on fine and coarse grids are both based on CIP-Pl. 



Level 


2 


3 


4 


5 


DOFs 


16641 


66049 


263169 


1050625 


iter (Case 1) 


52 


64 


66 


63 


iter (Case 2) 


27 


24 


22 


21 



shapes and amplitudes as the exact solutions. Table [3] presents the comparing of shifted Laplace 
operator with FEM-Pl and the original Helmholtz operator with CIP-Pl. For the case k = 100, we 
can see that the second approach has the minimum iteration counts. The iteration counts shown 
in Table H] examine the performance of PGMRES based on Algorithm 13.11 with different steps of 
GMRES smoothing. The CIP-Pl and CIP-P2 are only used for approximations on coarse grids. 
From Tables [TJ [3] and HI it suggests that Algorithm 13.11 is also very efficient as Algorithm 12.11 and 
the second approach in Table [3] which uses CIP-FEM on both fine and coarse grids, and the growth 
in GMRES smoothing steps does not always reduce iteration counts obviously. This supports the 
similar phenomena shown in the right graph of Figure [5l Hence, in the following we will only use 
one step (m2 = 1) of GMRES smoothing in Algorithm 13. 1[ 

Table 4: Iteration counts of PGMRES based on Algorithm [XT] for discrete system from FEM-Pl and FEM-P2 with 
different steps of post GMRES smoothing (k — 100). 



Level 


3 


4 


5 


DOFs 


66049 


263169 


1050625 


iter (PI, m2 = 1) 


24 


23 


22 


iter (PI, 7712 = 10) 


21 


23 


23 


iter (P2, 7712 = 1) 


27 


21 


19 


iter (P2, 7712 = 10) 


32 


21 


18 



Table 5: Iteration counts of PGMRES based on Algorithm [3TT] for discrete system from FEM-Pl and FEM-P2 for 
the cases k. = 50, 200, 360 with the coarsest grid condition nho/p ~ 2. 



K = 50 


Level 


3 


4 


5 


DOFs 


16641 


66049 


263169 


iter (PI) 


15 


15 


15 


iter (P2) 


23 


14 


13 


«: = 2QQ 


Level 


3 


4 


5 


DOFs 


263169 


1050625 


4198401 


iter (PI) 


49 


57 


55 


iter (P2) 


44 


41 


36 


«: = 36Q 


Level 


2 


3 


4 


DOFs 


263169 


1050625 


4198401 


iter (PI) 


111 


115 


112 


iter (P2) 


41 


44 


39 



Tables [S] and E] show the iteration counts of PGMRES based on Algorithm 13. II with CIP-FEM 
used for correction only on coarse grids for discrete system from FEM-Pl and FEM-P2. For fixed 
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Table 6: Iteration counts of PGMRES based on Algorithm EH] for discrete system from FEM-Pl and FEM-P2 
for the cases k — 400, 500, 600. The coarsest grid condition is chosen with mesh size ho ~ 0.00276 such that 
Kho/p ~ 1.1 ~ 1.7. 



FEM-Pl 


Level 


2 


3 


FEM-P2 


Level 


2 


3 


DOFs 


1050625 


4198401 


DOFs 


1050625 


4198401 


iter (k = 400) 


30 


26 


iter (k = 400) 


21 


14 


iter (k = 500) 


71 


50 


iter (k = 500) 


29 


25 


iter (k = 600) 


156 


141 


iter (k = 600) 


43 


47 



K, one can observe that the iteration counts are robust on different levels. Due to the same coarsest 
grid used for k = 400, 500, 600, the growth of iteration counts is a little faster than linear when 
K > 400. But when the FEM-P2 is applied in particular, the growth of iteration counts with 
increasing wave number is more stable than FEM-Pl. 

Example 5.2. In a unit square domain $7 with center (0,0), we consider the Helmholtz problem 
()l.ip - ()1.2p with discontinuous wave number, which is defined by 

f Ki, if (xi,X2) G (-0.5,0) X (0,0.5)U(0,0.5) x (-0.5,0), 

K = < 

I K2, elsewhere, 

where K2 = qKi,q > 1. We set the Robin boundary condition ()1.2p with g = and the external 
force f{x) to be a narrow Gaussian point source (cf. [7J) located at (ri,r2) = (-0.25,-0.25): 



Imaginary part of discrete FEM-P2 Solution 




Imaginary pari ol discrete FEM-P2 Sdjtion 





Figure 7: Surface plot of imaginary part of discrete FEM-P2 solutions for K2 
(right) on the grid with mesh condition K2h/p ~ 0.8. 



300, g = 3 (left) and K2 = 300, <j = 10 



In this example, we test Algorithm 13.11 for the Helmholtz problem with discontinuous wave 
number. Due to the fact that K2 = max{«;i,K2}, the coarsest mesh condition is according to the 
choice of K2ho/p. Figure [7| displays the surface plot of imaginary part of discrete FEM-P2 solutions 
for K2 = 300 with q = 3 and g = 10 on the grid with mesh condition K2h/p ~ 0.8. The iteration 
counts of PGMRES based on Algorithm 13.1! for different K2 and q are listed in Table [71 For fixed 
K2 and polynomial order, we note that the iteration counts are robust on different levels. Similar 
to the first example, the PGMRES iteration for FEM-P2 is more stable than FEM-Pl. 
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Table 7: Iteration counts of PGMRES based on Algorithm [5TT] for discrete system from FEM-Pl and FEM-P2 for 
the cases K2 = 180 (q = 3, 10), K2 = 300 (q = 3, 10). 



K2 = 180 


Level 


3 4 


DOFs 


263169 1050625 


iter (Pl,g = 3) 


26 27 


iter (Pl,g = 10) 


28 29 


iter (P2,g = 3) 


16 15 


iter {P2,q = 10) 


17 15 


K2 = 300 


Level 


3 4 


DOFs 


1050625 4198401 


iter (Pl,g = 3) 


30 30 


iter {Pl,q = 10) 


32 33 


iter (P2,g = 3) 


16 15 


iter (P2,g = 10) 


16 15 
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